We want to solve the linear system given by \begin{equation} Ax = b, \end{equation} where $A$ is a non-singular $n\times n$ matrix.
You have already solved this problem using Gaussian Elimination (and its partial pivoted version) which has an assymptotic cost $\mathcal{O}(n^3)$.
In this homework you will try to solve the system in lower complexity using an iterative method. In particular, you will implement the Jacobi and Gauss-Seidel iterations and you will study its limitations.
Both Jacobi and Gauss-Iterations can be seen as fixed point methods, used by a matrix splitting.
Given a square matrix $A$ you will define a splitting in two matrices $N$ and $M$ such that $A = N+M$. If you suppose that $N$ is invertible you can write the system \begin{equation} Ax = b, \end{equation} as \begin{equation} Nx = b - Mx, \end{equation} and you can define the fixed point iteration as \begin{equation} Nx^{n+1} = b - Mx^{n}, \end{equation} or equivalently \begin{equation} x^{n+1} = N^{-1}\left ( b - Mx^{n} \right), \end{equation} in which the inverse of $N$ is never computed explicitily.
In this case you can show that the convergence speed is can be determined from the spectral radius of the iteration matrix: \begin{equation} T = N^{-1}M. \end{equation}
If we define \begin{equation} f(x) := N^{-1}\left ( b - Mx^{n} \right), \end{equation} then we are solving the fixed point problem \begin{equation} x = f(x). \end{equation}
We know that in this case, the fixed point iteration given by $x^{n+1} = f(x^n)$ converges if and only if there exist $\kappa <1$ such that \begin{equation} \| f(x) - f(y) \| \leq \kappa \|x - y\|, \qquad \forall x,y \in R^n, \end{equation} where we are using the Euclidean metric given by \begin{equation} \| x \| = \sqrt{\sum_{i = 1}^n |x_i|^2 }. \end{equation}
In such case, the error is given by \begin{equation} \| x^{n+1} - x\| \leq \frac{\kappa^{n+1}}{1 - \kappa} \| x^{0} - x\|, \end{equation} in other words, the convergence and its ratio ar given by $\kappa$.
Using the expression for $f$ we have that \begin{align} \| f(x) - f(y) \| & = \| N^{-1}(b - Mx) - N^{-1}(b - My) \|, \\ & = \| N^{-1} M(x -y) \|, \\ & \leq \|N^{-1} M \| \|x - y \|; \end{align} or \begin{equation} \kappa = \|N^{-1} M \| = \|T \|. \end{equation}
When we are using the Euclidean norm, we have that \begin{equation} \|T \| = \max_{z \neq 0} \frac{\|Tz \|}{\| z\|} = \rho(T). \end{equation} where the spectral radius $\rho(T)$ is given by \begin{equation} \rho(T) = \max_{i = 1,..,n} |\lambda_i|, \qquad \mbox{ where } \lambda_i \mbox{ are the eigenvalues of } T \end{equation}
The Jacobi iteration correspond to a fixed point method, in which the Matrix splitting is given by \begin{equation} N = D, \qquad M= L+U, \end{equation} where $D$ is the diagonal of $A$ and $L$ and $U$ are the (strictly) upper and lower triangular parts of $A$.
In order to extract the diagonal, upper triangular and lower triangular matrices from $A$ we will use the built-in functions in Julia.
In [8]:
# size of the matrices
m = 10
# generate a random matrix
A = rand(m,m) + m*eye(m)
# get the digonal part
D = diagm(diag(A),0)
# to get the upper triangular part of the matrix
U = triu(A,1)
#and the lower part
L = tril(A,-1)
# we check that this is indeed a matrix splitting
println("Error in the splitting = ",norm(A -(L+D+U)))
For debugging and testing your functions you may find the macro @show usefull, if you type @show at the beginning of a line, it will force julia to show in the console the command, and what it is doing.
Q1.a Implement the Jacobi iteration with input:
Your function will have the following optional parameters
The ouput of your function is the final approximation $x$ of your vector, of in the case that history is true, it will output a matrix with the all the intermediate approximations of size $ n \times \mbox{number of iterations to converge}$.
Your function will raise an error if it didn't converge in the $Nmax$ iterations, and if the dimension of $A$ and $b$ are not consistent.
Hint:
In [ ]:
function JacobiIt(A,b; Nmax = 30, ϵ=1e-6, history = false)
end
The Gauss-Seidel iteration is analogous to the Jacobi iteration; however, it uses a different splitting, namely \begin{equation} N = D+U, \qquad M= L, \end{equation}
Q1.a Implement the Gauss-Seidel iteration with input:
Your function will have the following optional parameters
The ouput of your function is the final approximation $x$ of your vector, of in the case that history is true, it will output a matrix with the all the intermediate approximations of size $ n \times \mbox{number of iterations to converge}$
Your function will raise an error if it didn't converge in the $Nmax$ iterations, and if the dimension of $A$ and $b$ are not consistent.
Hint:
In [ ]:
function GaussSeidelIt(A,b; Nmax = 30, ϵ=1e-6, history = false)
end
In this question you will need to compute the spectral radius of matrices, for a matrix $A$, its spectral radius $\rho (A)$ can be computed by
In [2]:
# generating random matrix of size m \times m
m = 100
A = rand(m,m) + m*eye(m)
# computing the eigenvalues and eigenvectors of A
(λ, X) = eig(A)
# obtaining the one with the maximum absolute value
rhoA = maximum(abs(λ))
Out[2]:
Use your algorithms to solve the following randomly generated system
In [ ]:
m = 100
A = rand(m,m) + m*eye(m)
b = rand(m,1)
@time XJac = JacobiIt(A,b, history = true)
println("The residual is : ", norm(A*XJac[:,end]- b)/norm(b) )
println("number of iterations is : ",size(XJac,2))
@time XGS = GaussSeidelIt(A,b, history = true)
println("The residual is : ", norm(A*XGS[:,end]- b)/norm(b) )
println("number of iterations is : ",size(XGS,2))
Q3.a How can you explain the one converges faster than the other? Write a small script that computes the spectral radius of the iteration matrix for both algorithms and use it on your explanation.
In [ ]:
# write your script in here
Answer
Now you will use your algorithms to solve a slighly different system
In [ ]:
m = 100;
A = rand(m,m) + sqrt(m)*eye(m);
b = rand(m,1);
With the Jacobi iteration,
In [ ]:
@time XJac = JacobiIt(A,b, history = true)
println("number of iterations is : ",size(XJac,2))
and with the Gauss Seidel iteration
In [ ]:
@time XGS = GaussSeidelIt(A,b, history = true)
println("number of iterations is : ",size(XGS,2))
Q3.b How can you explain the one converges and the other just fails? Write a small script that computes the spectral radius of the iteration matrix for both algorithms and use it on your explanation. Hint:
In [ ]:
# write your script in here
Answer:
Q4.a What is the complexity of one iteration of the Gauss-Seidel method? what about complexity of the one iteration of the Jacobi? What would be the complexity the matrix was sparse? Can you easily parallelize the Jacobi iteration? what about Gauss Seidel?
Answer:
Q4.b What would be the condition on the number of iterations such that the Jacobi and Gauss-Seidel iteration would have a better complexity than Gaussian Elimination?
Answer:
Q4.c Run a benchmark with the randomly generated systems that you use to test your functions. Can you achieve a complexity $\mathcal{O}(n^2)$?
In [ ]:
# write your whole script here